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Abstract 

In this paper, we provide a matrix-analytic solution for randomized load balanc- 
ing models (also known as supermarket models) with phase-type (PH) service times. 
Generalizing the service times to the phase-type distribution makes the analysis of 
the supermarket models more difficult and challenging than that of the exponential 
service time case which has been extensively discussed in the literature. We first de- 
scribe the supermarket model as a system of differential vector equations, and provide 
a doubly exponential solution to the fixed point of the system of differential vector 
equations. Then we analyze the exponential convergence of the current location of 
the supermarket model to its fixed point. Finally, we present numerical examples to 
illustrate our approach and show its effectiveness in analyzing the randomized load 
balancing schemes with non-exponential service requirements. 



1 Introduction 

In the past few years, a number of companies (e.g., Amazon, Google, ,...etc) are offering 
the cloud computing service to enterprises. Furthermore, many content publishers and 
apphcation service providers are increasingly using Data Centers to host their services. 
This emerging computing paradigm allows service providers and enterprises to concentrate 
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on developing and providing their own services/goods without worrying about computing 
system maintenance or upgrade, and thereby significantly reduce their operating cost. 
For companies that offer cloud computing service in their data centers, they can take 
advantage of the variation of computing workloads from these customers and achieve the 
computational multiplexing gain. One of the important technical challenges that they 
have to address is how to utilize these computing resources in the data center efficiently 
since many of these servers can be virtualized. There is a growing interest to examine 
simple and robust load balancing strategies to efficiently utilize the computing resource of 
the server farms. 

Distributed load balancing strategies, in which individual job (or customer) decisions 
are based on information on a limited number of other processors, have been studied ana- 
lytically by Eager, Lazokwska and Zahorjan [HEKG] and through trace-driven simulations 
by Zhou [26j. Further, randomized load balancing is a simple and effective mechanism 
to fairly utilize computing resources, and also can deliver surprisingly good performance 
measures such as reducing collisions, waiting times, backlogs,... etc. In a supermarket 
model, each arriving job randomly picks a small subset of servers and examines their 
instantaneous workload, and the job is routed to the least loaded server. When a job 
is committed to a server, jockeying is not allowed and each server uses the first-come- 
first-service (FCFS) discipline to process all jobs, e.g., see Mitzenmacher [11^ I12j. For the 
supermarket models, most of recent research applied density dependent jump Markov pro- 
cesses to deal with the simple case with Poisson arrival processes and exponential service 
times, and illustrated that there exists a fixed point which decreases doubly exponentially. 
Readers may refer to, such as, a simple supermarket model by [T| 124 ^ [TT | [T2j : simple varia- 
tions bydSmamilTlEailElES]; bad information by [201 El [ISl [IS]; fast Jackson network 
by Martin and Suhov [10l [9l[2T]: and general service times by Bramson, Lu and Prabhakar 
[2]. When the arrival processes or the service times are more general, the available results 
of the supermarket models are few up to now. The purpose of this paper is to provide a 
novel approach for studying a supermarket model with PH service times, and show that 
the fixed point decreases doubly exponentially. 

The remainder of this paper is organized as follows. In the next section, we describe the 
supermarket model with the PH service times as a system of differential vector equations 
based on the density dependent jump Markov processes. In Section 3, we set up a system 
of nonlinear equations satisfied by the fixed point, provide a doubly exponential solution to 
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the system of nonlinear equations, and compute the expected sojourn time of any arriving 
customer. In Section 4, we study the exponential convergence of the current location 
of the supermarket model to its fixed point. In Section 5, numerical examples illustrate 
that our approach is effective in analyzing the supermarket models from non-exponential 
service time requirements. Some concluding remarks are given in Section 6. 

2 Supermarket Model 

In this section, we describe a supermarket model with the PH service times as a system 
of differential vector equations based on the density dependent jump Markov processes. 

Let us formally describe the supermarket model, which is abstracted as a multi-server 
multi-queue stochastic system. Customers arrive at a queueing system of n > 1 servers as 
a Poisson process with arrival rate nX for A > 0. The service times of these customers are 
of phase type with irreducible representation (a, T) of order m. Each arriving customer 
chooses d > 1 servers independently and uniformly at random from these n servers, and 
waits for service at the server which currently contains the fewest number of customers. 
If there is a tie, servers with the fewest number of customers will be chosen randomly. All 
customers in every server will be served in the FCFS manner. Please see Figure 1 for an 
illustration. 

Poisson Arrival ^--Y^ [ | :| Q 

nx — znmo 

each customer 
probes d servers 

mm ® 

Figure 1: The supermarket model: each customer can probe the loading of d servers 

For the supermarket models, the PH distribution allows us to model more realistic sys- 
tems and understand their performance implication under the randomized load balancing 
strategy. As indicated in [7], the process lifetime of many parallel jobs, in particular, 
jobs to data centers, tends to be non-exponential. For the PH service time distribution, 
we use the following irreducible representation: (a,T) of order m, the row vector a is a 
probability vector whose jth entry is the probability that a service begins in phase j for 
1 < j < m; T is an m X m matrix whose {i,jY^ entry is denoted by tij with tj^j < 
for 1 < i < m, and tij > for 1 < i,j < m and i ^ j. Let = —Te ^ 0, where e is 
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a column vector of ones with a suitable dimension in the context. The expected service 
time is given by 1/^ = —aT~^e. Unless we state otherwise, we assume that all random 
variables defined above are independent, and that the system is operating in the stable 
region p = X/ fi < 1. 

We define n^*^ (t) as the number of queues with at least k customers and the service 
time in phase i at time t > 0. Clearly, < n^*"* (t) < n for k > and 1 < i < m. Let 



n 



n 



1, 



and k > 1 



n 



which is the fraction of queues with at least k customers and the service time in phase i 
at time t > 0. We write 

x^J^) it) = [xt''^ (t) , xlt^^) (t) , . . . , 4^'-) it)) , > 1, 



The state of the supermarket model may be described by the vector Xn {t) for t > 0. 
Since the arrival process to the queueing system is Poisson and the service times of each 
server are of phase type, the stochastic process {Xn (t) ,t > 0} describing the state of the 
supermarket model is a Markov process whose state space is given by 



Let 



and k > 1 



and ng^^ is a vector of nonnegative integers for k > 1}. 



t) = E \xt^^ it) 



Clearly, sq {n,t) = 1. We write 



k (n, t) = (4'^ (n, t) , (n, t) , . . . , s^^^ (n, i)) , k> 1. 



As shown in Martin and Suhov [TO] and Luczak and McDiarmid [8], the Markov process 
{Xn (t) ,t > 0} is asymptotically deterministic as n — )• oo. Thus the limits lim„_s>oo E Xn^ (t) 
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and lim„^oo 
we write 

for A; > 1 



E 



X 



always exist by means of the law of large numbers. Based on this, 

{t) = lim So (n,t) = 1, 



sf it) = lim s« 
Sk{t) = {s^k'' (t),4'^ it) 



and 



S{t) = {So{t),S, {t),S2{t),...). 



Let X {t) = lim„_^oo -^n Then it is easy to see from the Poisson arrivals and the PH 
service times that {X (t) ,t> 0} is also a Markov process whose state space is given by 

^ = { {9^'\9^'\9^'\ ...): 9^'^ = l,9^'-''> > 9^'^ > o} ■ 

If the initial distribution of the Markov process {Xn (t) ,t > 0} approaches the Dirac 
dclta-mcasurc concentrated at a point g £ ^1, then its steady-state distribution is con- 
centrated in the limit on the trajectory Sg = : t > 0}. This indicates a law of 
large numbers for the time evolution of the fraction of queues of different lengths. Fur- 
thermore, the Markov process {Xn (t) ,t > 0} converges weakly to the fraction vector 
S (t) = {So {t) ,Si{t) ,S2{t) ,...), or for a sufficiently small £ > 0, 

lim P{\\Xnit)-S{t)\\>e} = 0, 

n— ^-oo 

where ||a|| is the Loo-norm of vector a. 

In what follows we provide a system of differential vector equations in order to de- 
termine fraction vector S{t). To that end, we introduce the Hadamard Product of two 
matrices A = {aij) and B = as follows: 



Specifically, for A; > 2, we have 



AQB = {aijbij). 



A®'' = Aq Aq ■ ■ ■ Q A. 



k matrix A 

To determine the fraction vector S (t) , we need to set up a system of differential vector 
equations satisfied by S {t) by means of the density dependent jump Markov process. To 
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that end, we provide a concrete example for A; > 2 to indicate how to derive the the system 
of differential vector equations. 

Consider the supermarket model with n servers, and determine the expected change 
in the number of queues with at least k customers over a small time period of length dt. 
The probability vector that during this time period, any arriving customer joins a queue 
of size — 1 is given by 



n 



XS^^,in,t)-XSfin,t) 



dt. 



Similarly, the probability vector that a customer leaves a server queued by k customers is 
given by 

n [Sk (n, t)T + Sk+i (n, t) r°a] dt. 



Therefore we can obtain 



dE [nk (n, t)] =n XSf^, (n, t) - XS'^" (n, t) 



dt 



+ n [Sk {n, t)T + Sk+i (n, t) T^a] dt, 

which leads to 

^^^^1^ = XS®^^ (n, t) - XSf (n, t) + (n, t)T + Sk+i (n, t) T'^a. 

Taking n — >■ oo in the both sides of Equation (??), we have 

= A^'i (i) - ^Sf {t) + Su {t) T + Su^, it) T^a. 

Using a similar analysis to Equation (??), we can obtain a system of differential vector 
equations for the fraction vector S (t) = {Sq (t) , (t) , 52 (t) , . . .) as follows: 



dt 



So (t) = 1, 
So{t) = -XS^{t) + Si{t)T^, 



-Si {t) = XaS^ {t) - XSf^ (t) + Si {t) T + S2 (t) T^a, 



and for k >2, 



dt 



Sk {t) = A5®_^ (t) - XSf {t) + Sk (t) T + Sk+i (t) T^a. 



(1) 
(2) 

(3) 
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Remark 1 Mitzenmacher 12^ provided an heuristical and interesting method to estab- 
lish such systems of differential equations, but they lack a rigorous mathematical meaning 
for understanding the stochastic process {Xn {t),t> 0} in which Xn {t) = {X^n^ (t) , Xi^^ (t) , 
Xn^ (t) , . . .) and Xn^ (t) = (t) /n for k > 0. This section, following Martin and Suhov 
111 Of and Luczak and McDiarmid 18], gives some necessary mathematical analysis for the 
stochastic process {Xn (t) ,t > 0} and the system of differential vector equations (OP, 
and 



3 A Matrix- Analytic Solution 

In this section, we provide a doubly exponential solution to the fixed point of the system 
of differential vector equations ([T|) , ([2|) and ^ . 

A row vector vr = (ttq, tti, 7r2, . . .) is called a fixed point of the fraction vector S (t) if 
limf_j.+oo S (t) = vr. In this case, it is easy to see that 



lim 

i— i-+oo 



0. 



Therefore, as t ^ +oo the system of differential vector equations ([I]), ([2]) and ([3]) can be 
simplified as 

- A4 + ^iT" = 0, (4) 
Aa4 - Avrf^ + ttiT + iiiT^a = 0, (5) 



and for k > 2, 



Avr®f 1 - Avrf + vr^T + ^fc+iT^Q = 0. (6) 



In general, it is more difficult and challenging to express the fixed point of the super- 
market models with more general arrival processes or service times, because the systems 
of nonlinear equations are more complicated for computation. Fortunately, we can derive 
a closed- form expression for the fixed point vr = (ttq, vri, 7r2, ...) for the supermarket model 
with PH service times by means of a novel matrix-analytic approach given as follows. 

Noting that Sq (t) = 1 for all t > 0, it is easy to see that ttq = 1. It follows from 
Equation dH that 

ttiTO = A. (7) 
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To solve Equation ^ , we denote by uj the stationary probability vector of the irreducible 
Markov chain T + T^a. Obviously, we have 

-uT^ = A. (8) 

Thus, we obtain vri = ^oj = p ■ oj. Based on the fact that ttq = 1 and tti = p • it follows 
from Equation ([5|) that 

\a - a/ ■ uj'^'^ + p ■ ujT + 712 r°a = 0, 

which leads to 

A - a/ • uj^'^e + p • ioTe + TrsT^ = 0. 
Note that ojTe = —p, we obtain 

vraTO = Xp'^uj^^'^e. 
Let = oj'^'^e. Then it is easy to see that 9 G (0, 1), and 

vraT^ = XOp"^. 
Using a similar analysis to Equation ([8]), we have 

^^ = ^u^ = epd+K^, (9) 

Based on vri = p ■ ui and 7r2 = Op'^^^ ■ u, it follows from Equation ([6]) that for k = 2, 
which leads to 

A^/ - A^'^+V'^'^'^ + ^Z'*'^ • wTe + ^3^° = 0, 

thus we obtain 

ttsTO = A0'^+i/+''. 
Using a similar analysis on Equation dH), we have 

= ^^^-^ a; = ^'i+iZ+rf+i . a;. (10) 



Based on Equations ([9]) and pOj) . we may infer that there is a structured expression 
TT^ = 0'^'° i-'^+ip'^'' ^+"''' ^^ *"°'+i-6t;, for > 1. To that end, the following theorem 

states this important result. 
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Theorem 1 The fixed point vr = (ttq, tti, 7r2, . . .) is unique and is given by 



Vro = 1, TTl = p-u 

and for k > 2, 

gd'-— 2+(i'=-3H hi ^d'^-i+d'-'-^H hi 



vTfc = 6'' +^ ^+-+ip'i^ ^+<i- . (11) 



or 



Proof: By induction, one can easily derive the above result. 

It is clear that Equation (jlip is correct for the cases with / = 2, 3 according to Equations 
([9]) and (jlOp . Now, we assume that Equation (jlip is correct for the cases with I = k. Then 
it follows from Equation ^ that for / = /c + 1, we have 

^Qd''-^+df'-''^+-+dpd'''^+d'''^+ - +d . ^Qd _ y^Qd^-^+d''-^ + ---+d pd^+d^-^ + ---+d . ^Qd 

which leads to 



thus we obtain 



By a similar analysis to ([8]), we have 

y^Qd'^-^+d''''^^ ^d+l pd'^+d''^^^ Vd 

_ £ifc+(ifc-i_| |_rf4.i 

This completes the proof. I 
Now, we compute the expected sojourn time that a tagged arriving customer spends 
in the supermarket model. For the PH service times, a tagged arriving customer is the kth. 
customer in the corresponding queue with probability vector vr^^^ — vr^''. When A; > 1, 
the head customer in the queue has been served, and so its service time is residual and is 
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denoted as Xr. Let X be of phase type with irreducible representation (a,r). Then Xr 
is of phase type with irreducible representation {u^T). Clearly, we have 



E[X] = a{-T)-^e, E[Xr] 



-T)-^e. 



Thus it is easy to see that the expected sojourn time of the tagged arriving customer is 
given by 

oo 

E [T,] = (4 - '^e) E[X] + Y^ (vrf - J e {E [Xr] + kE [X]} 



k=l 



7Tfe{E[Xn]-E[X]} + E[X] 



1 + 



k=l 



p'^e {oj - a) {-ry^ e + a i-T)-^ 



k=l 



When the arrival process and the service time distribution are Poisson and exponential, 
respectively, it is clear that a = uj = 9 = 1 and a {—T)~^ e = thus we have 

y p d-1 



Em 



^ oo 



k=0 



which is the same as Corollary 3.8 in Mitzenmacher [12]. 

In what follows we consider an interesting problem: how many moments of the service 
time distribution are needed to obtain a better accuracy for computing the fixed point or 
the expected sojourn time. It is well-known from the theory of probability distributions 
that the first three moments is basic for analyzing such an accuracy, and we can construct 
a PH distribution of order 2 by using the first three moments. Telek and Heindl [22] 
provided a fitting procedure for matching a PH distribution of order 2 from the first three 
moments exactly. It is necessary to list the fitting procedure as follows: 

For a nonnegative random variable X, let m„ = n > 1. We take a PH 

distribution of order 2 with the canonical representation (a, T), where a = {t],1 — rj) and 




< r/ < 1 and < < ^2- Note that the three unknown parameters rj, and ^2 can be 
obtained from the first three moments mi, m2 and of an arbitrary general distribution. 

In Table 1, = m2/m\ — 1 is the squared coefficient of variation. If the moments 
do not satisfy these conditions in Table 1, then we may analyze the following four cases: 
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Tabic 1: Specific Bounds of the First Three Moments 



Moment 


Condition 


Bounds 


mi 




< mi < oo 


m2 




1.5ml — "^2 




0.5 < 4 < 1 


3/?/f - 1 + v/2 (1 - 4)^) < < G"'?4 




i<4 


|mf (1 + 4)^ < ma < (X) 



ma 



(a.l) if m2 < 1.5mf , then we take m2 = 1.5mf; 

(a.2) if 0.5 < 4 ^ 1' and ma < 3m^ (^34 - 1 + (l - 4) ^) ' then we take 

3mf (34-1 + ^/2(1-4)^); 

(a.3) if 0.5 < 4 — "^3 > 6m;|4) then we take ma = 6mf4; ^.nd 

(a.4) if 1 < 4) and ma < |mf (l + 4)^' then we take ma = |mf (l + 4)^- 

Let c = SmI — 2mima, d = 2mf — m2, 6 = 3mim2 — ma and a = — 6cd. If the 

moments respectively satisfy their specific bounds shown in Table 1 or the exceptive four 

cases, then three unknown parameters rj, ^1 and ^2 can be computed in the following three 

cases. 

(1) If 0, then 



—b + 6mld + ^/a h — JH b + ^/a 
V = J— — 7= ' ^1 = ' « = ■ 



b + ^/a 



(2) If c < 0, then 



b — 6mld + ^/a b + ^/a b—^fa 
V = — ]= > ?i = > ?2 = • 



(3) If c = 0, then 



-b + ^/a 



r? = 0, 6 > 0, 6 = — • 
mi 



Prom the above discussion, we can always construct a PH distribution of order 2 to 
approximate an arbitrary general distribution with the same first three moments. In fact, 
such an approximation achieves a better accuracy in computation. 

For the PH distribution of order 2, we have 




r) 1 — r) 



H 



^2V -^2V 
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which leads to 



UJ = 



) 



and 



e 



Thus, the PH distribution of order 2 can effectively approximates an arbitrary general 
service time distribution in the supermarket model under the same first three moments, 
and specifically, all the computations are very simple to implement. 

Remark 2 Bramson, Lu and Prabhakar JS] provided a modularized program based on 
ansatz for treating the supermarket model with a general service time distribution. They 
organized a functional equation tt = F {G (tt)) for analyzing the stationary probability vec- 
tor IT in terms of insensitivity and generalized Fibonacci sequences, although the operators 
F and G are not easy to be given for this supermarket model. This paper studies the super- 
market model with a PH service time distribution, provides the doubly exponential solution 
to the fixed point, and is specifically related to the phase type environment by means of 
the crucial factor 9 = a;®'^e. Note that the PH distributions are dense in the set of all 
nonnegative random variables, this paper can numerically provide necessary understanding 
for the role played by the general service time distribution in performance analysis of the 
supermarket model by means of the PH approximation of order 2. 

4 Exponential convergence to the fixed point 

In this section, we study the exponential convergence of the current location S (t) of the 
supermarket model to its fixed point vr. 

For the supermarket model, the initial point S (0) can affect the current location S (t) 
for each t > 0, since the service process in the supermarket model is under a unified 
structure. To that end, we provide some notation for comparison of two vectors. Let 
a = (ai, 02, 03, . . .) and 6 = (61, b2, 63, • • •). We write a -< 6 if < 6^ for some k > 1 
and ai < hi for I ^ k,l > 1; and a ^ b if Ok < bk for all k > 1. Now, we can obtain the 
following useful proposition whose proof is clear from a sample path analysis and thus is 
omitted here. 



Proposition 1 If S (0) ^5(0), then S (t) ^S{t). 
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Based on Proposition [H the following theorem shows that the fixed point vr is an upper 
bound of the current location S (t) for all t > 0. 

Theorem 2 For the supermarket model, if there exists some k such that Sk (0) = 0, then 
the sequence {S^ (t)} has an upper bound sequence which decreases doubly exponentially 
for all t > 0, that is, S (t) ^ ir for all t >0. 

Proof: Let Sk (0) = vr^. for A; > 1. Then for each k > 1, Sk{t) = Sk (0) = vrfc for all 
t > 0, since S (0) = (0) , ^2 (0) , 52 (0) , . . is a fixed point in the supermarket model. 
If Sk (0) = for some k, then Sk (0) -< Sk (0) and Sj (0) ^ Sj (0) for j ^ k,j >1, thus 
S (0) ^ S (0). It is easy to see from Proposition [1] that Sk (t) ^ Sk (t) = iTk for all > 1 
and t > 0. Thus we obtain that for all > 1 and t > 

Sk (t) <d'i-i p d-1 . cj. 

This completes the proof. I 
To show the exponential convergence, we define a Lyapunov function $ (t) as 

oo 

$ (t) = ^ Wk [iTk - Sk (t)] e 
k=l 

in terms of the fact that 5^ (t) ^ vr^ for k > 1 and ttq = So{t) = 1, where {wk} is a 
positive scalar sequence with Wk+i > li'fc > if i = 1 for k > 2. 

The following theorem measures the distance $ (t) of the current location S (t) for 
t > to the fixed point vr, and illustrates that this distance between the fixed point and 
the current location is very close to zero with exponential convergence. This shows that 
from a suitable starting point, the supermarket model can be quickly close to the fixed 
point. 

Theorem 3 For t > 0, $ (t) < coe^*^*, where cq and 6 are two positive constants. In this 
case, the potential function ^ (t) is exponentially convergent. 

Proof: Note that 

oo 

= J2wk [TTk- Sk{t)]e, 
k=l 

we have 

^^Ht) = -±n^kpkit)e. 

k=l 
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It follows from Equations ([T]) to ([3|) that 

(t) = - wi[\S^ {t) a - XSf (t) + Si (t) T + S2 (t) r°a]e 



J2 Wk[>^Sf^i (t) - ^Sf (t) + Sk (t) T + Sk+i (t) rOQ]e. 



k=l 



By means of Sq (t) = 1 and Te = —T^, we can obtain 

^^^(t) = -wi[X- XSf (t) e - Si (t) ro + 52 it) T^] 



M>^Sf^i (t) e - XS'i'' (t) e - Sk (t) + 5,+i (t) 



k=2 

We take some nonnegative constants {t) and dj. {t) for k>l such that 



X = fi{t) Si {t)T\ 

for k>l 

XSf (t) e = Ck (t) [TTk - Sk it)] e 

and 

5fe(t)rO = 4 (t) [TTk-Sk{t)]e. 
Then it follows from (|12|) that 

it) = - {[iw2 - Wl)] Cl (t) + Wi [fi it) - 1] dl (t)} ■ [TTi - Si (t)] e 



^ [{wk+i - Wk) Ck (t) + {wk-i - Wk) dk (t)] • [vTfc - Sk (t)] e. 



k=2 

For a constant 5 > 0, we take 

Wl = 1, 

[(i/;2 - Wl)] Cl (t) + Wl [fi (t) - 1] dl (t) > 6wi 

and 

{wk+i - Wk) Ck (t) + (wk-i - Wk) dk {t) > 6wk. 
In this case, it is easy to see that 

W2>1 + 7- 

Cl (t) 

and for k > 2 

6wk dl: (i) , N 
Wk^i >Wk-\ 7— H 7— [Wk - Wk-i) . 

Ck [t) Ck [t) 
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Thus we have 



which can leads to 



— $ (t) < -5 ^ Wk [TTfc - Sk (t)] e = -5$ {t) , 



fe=0 



$ it) < coe 



-St 



This completes the proof. 



5 Numerical examples 



In this section, we provide some numerical examples to illustrate that our approach is 

effective and efficient in the study of supermarket models with non-exponential service 
requirements, including Erlang service time distributions, hyper-exponential service time 
distributions and PH service time distributions. 

Example one (Erlang Distribution) We consider an m-order Erlang distribution with 
the irreducible PH representation (a, T), wherea = (1, 0, . . . , 0, 0) and 



T = 



-r) rj 

— T] T] 



\ 



—7] rj 



-1 1 








It is clear that 



T + T^a 



\ 



-rj ri 
\ V -1] J 

which leads to the stationary probability vector of the Markov chain T + T^a as follows: 

d 



; 1^ = coT"" = -L; p = - = ; 9 = mi-] =m^'^. 

\m m m mj m fx r] \m J 
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Thus we obtain 



V V 




d-1 



1 1 



1 1 



TTk=m 



m m 



m m 



1 1 



1 1 



= m 



m m 



m m 



Let A = 1. If p = ^ < 1, then this supermarket model is stable. In the stable case, 
r] > m. We may consider the following simple cases: 

(a) If m = 2 and d = 2, then vr^ = 2^'"'' r]^'^" . 

(b) If m = 3 and d = 2, then iTk = S^'^'r/i-^'. 

Based on the two simple examples with A = 1 and d = 2, we need to illustrate how 
the fixed point depends on the stage number m and the exponential service rate r]. To 
that end, we write Tr^ {m,r]). It is easy to see that for a given pair {k,r]) for r] > m and 
= 1, 2, . . . , we have 



On the other hand, for a given pair {k, m) for m. A; = 1, 2, . . . , we can see that vr^ (m, rj) is 
a decreasing function of tj. 

Let us consider the average response time of the supermarket model with an m— stage 
Erlang distribution. We first consider a parallel system with n = 100 servers and the 
service time distribution is exponential. We normalize the average service time to unity 
and vary the arrival rate A. For the m— stage Erlang distribution, the bigger the number m 
is, the bigger its variance is. Table [2] illustrates the average response time under different 
probe size d. One can observe that there is a dramatic improvement (or reduction) in the 
average response time when increasing the probe size d. 

We further analyze the cases that the service time is either distributed according to 
2-stage Erlang or 3-stage Erlang distribution. Similarly, we normalized the total average 
service time as unity and we vary the arrival rate A. Tables [3] and H] illustrate the average 
response time under different probe size d. One can observe that 

• Simple probing size d can significantly improve the performance by lowering the 
average response time. 

• When the service time has lower variance, the average response time is lower. 



TTk (1, 7?) < VTfc (2, r/) < • • • < VTfc (m, T/) < • • • . 
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Table 2: Average response time for exponential service time 



number of servers (n) 


probe size (d) 


arrival rate (A) 


response time (ii[T]) 


100 


2 


0.500000 


1.395977 


100 


2 


0.700000 


1.768194 


100 


2 


0.800000 


2.072020 


100 


2 


0.900000 


2.721852 


100 


3 


0.500000 


1.395320 


100 


3 


0.700000 


1.604113 


100 


3 


0.800000 


1.802933 


100 


3 


0.900000 


2.209601 


100 


5 


0.900000 


1.916280 



Table 3: Average response time for 2— stage Erlang service time 



number of servers (n) 


probe size (d) 


arrival rate (A) 


response time (B[T]) 


100 


2 


0.500000 


1.353783 


100 


2 


0.700000 


1.599851 


100 


2 


0.800000 


1.829199 


iOO 


2 


O.DOOOOO 


2.298170 


100 


3 


0.500000 


1.325610 


100 


3 


0.700000 


1.492651 


100 


3 


0.800000 


1.639987 


100 


3 


0.900000 


1.941196 


iOO 


T) 


O.IKKXKX) 


l.T:i9<s(iV 



Example two (Hyper-Exponential Distribution) We consider an m-order hyper-exponential 

m 

distribution F (x) = 1 — CKfe 6xp {— r^fex}, or the probability density function f (x) = 



k=l 



afc77fc exp {— r^fcx}. It is clear that the hyper-exponential distribution is of phase type 

k=l 

with the irreducible representation (a, T), where a = {ai,a2, ■ ■ ■ , am), and 



-m 



-m 



which lead to 



T + r"a = 



-•ni (1 - ai) 



\ 



-rim ) 

??1"2 
-7/2 (1 - OL2) 

'r]mOL2 



-Vm (1 - am) J 
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Table 4: Average response time for 3— stage Erlang service time 



number of servers (n) 


probe size (d) 


arrival rate (A) 


response time (EIT]) 


100 


2 


0.500000 


1.322544 


100 


2 


0.700000 


1.539621 


100 


2 


0.800000 


1.739972 


100 


2 


0.900000 


2.148191 


100 


3 


0.500000 


1.298863 


100 


3 


0.700000 


1.452785 


100 


3 


0.800000 


1.581663 


100 


3 


0.900000 


1.834704 


100 


5 


0.900000 


1.678233 



In general, the system of equations uj (T + T^a) = and ue = 1 does not admit a simple 
analytic solution. For a convenient description, we only consider a simple one with m = 2. 
In this case, we obtain 

Vim (ai + "2) 



and 



UJ 

A 



air/2 



a2m 



air/2 + 027?! ai??2 + Ci2Vl 



airj2 + a2Vi 



A(qi7?2 + «2??l) 

mm (ai + "2) 



air]2 



aim + a2??i 



+ 



a2m 



aim + a2m 



ai^2 
ai??2 + a2m 

A (air/2 + a2??i 
mm (ai + 02) 



+ 



a2m 



aim + a2m 



aim 



a2m 



^air]2 + a2m ai??2 + a2?/i. 
Tables [5] and [6] indicate how the doubly exponential solution (vri to vrs) depends on the 
vectors r] = [rji^m) and a = (01,02), respectively. 



Table 5: The doubly exponential solution depends on rj 





r, = (3,3) 


77 = (3, 10) 


J7= (3,20) 




(0.1667, 0.1667) 


(0.1667, 0.0500) 


(0.1667, 0.0250) 


7r2 


(0.0093, 0.0093) 


(0.0050, 0.0015) 


(0.0047, 0.0007) 


■KZ 


(2.858e-05, 2.858e-05) 


(4.626e-06, 1.388e-06) 


(3.819e-06, 5.728e-07) 


7I"4 


(2.722e-10, 2.722e-10) 


(3.888e-12, 1.166e-12) 


(2.485e-12, 3.728e-13) 




(2.470e-20, 2.470e-20) 


(2.746e-24, 8.238e-25) 


(1.053e-24, 1.579e-25) 



Let us consider the average response time of the supermarket model with an m-stage 
hyper-exponential service time distribution. We consider a parallel system with n = 100 
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Table 6: The doubly exponential solution depends on a 





Q = (0.5, 0.5) 


a = (0.2, 0.8) 


a = (0.8, 0.2) 


TTl 


(0.1667, 0.1667) 


(0.0667, 0.0267) 


(0.2667, 0.0067) 


^2 


(0.0047, 0.0005) 


(0.0003, 0.0001) 


(0.0190, 0.0005) 


TTS 


(3.680e-06, 3.680e-07) 


(9.136e-09, 3.654e-09) 


(9.607e-05, 2.402e-06) 


7r4 


(2.280e-12, 2.280e-13) 


(6.454e-18, 2.582e-18) 


(2.463e-09, 6.157c-ll) 




(8.752e-25, 8.752e-26) 


(3.221e-36, 1.289e-36) 


(1.618e-18, 4.046C-20) 



servers and the probability density function of the service time of a customer is given by 

fix) = 0.5 X (2 X + 0.25 x (0.5 x e""-^^) + 0.25 x (e"^). 

Note that the total average service time is normalized to unity and we vary the arrival 
rate A. Table [7] illustrates the average response time under different probe size d. One can 
observe that there is a dramatic reduction in the average response time when increasing 
the probe size. Furthermore, when the service time has a higher variance (we here compare 
it with the exponential distribution or m— stage Erlang distribution), the average service 
time is much higher. This indicates that we improve the performance of the supermarket 
model, one has to increase the probe size d. 

Table 7: Average response time for 3— stage Hyper-exponential service time 



number of servers (n) 


probe size (d) 


arrival rate (A) 


response time (£[?"]) 


100 


2 


0.500000 


1.552282 


100 


2 


0.700000 


1.969132 


100 


2 


0.800000 


2.360255 


100 


2 


0.900000 


3.225117 


100 


3 


0.500000 


1.462128 


100 


3 


0.700000 


1.723764 


100 


3 


0.800000 


1.947548 


100 


3 


0.900000 


2.476718 


100 


5 


0.900000 


2.066462 



Example three (PH Distribution) We consider an m-order PH distribution with irre- 
ducible representation (q,T). For m = 2,d = 2,a = (1/2, 1/2) and 



r(i) 




.T(2) 




,T(3) 




Table [8] illustrates how the doubly exponential solution depends on the PH matrices T (1), 
T (2) and T(3), respectively. 
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Table 8: The doubly exponential solution depends on the PH matrices T{i) 





T(l) 


T(2) 


T(3) 


TTl 


(0.2045, 0.1591) 


(0.1410, 0.1026) 


(0.3125, 0.2500) 


7r2 


(0.0137, 0.0107) 


(0.0043, 0.0031) 


(0.0500, 0.0400) 


713 


(6.193e-05, 4.817c-05) 


(3.965C-06, 2.884e-06) 


(0.0013 , 0.0010) 


714 


(1.259e-09, 9.793e-10) 


(3.390e-12, 2.465e-12) 


(8.446e-07, 6.757e-07) 


""5 


(5.204e-19, 4.048e-19) 


(2.478e-24, 1.802e-24) 


(3.656e-13, 2.925e-13) 




Figure 2: E [T^Js of the PH and exponential distributions for T(l) and T(2), respectively 

To discuss how different caused by a non-exponential distribution versus an exponen- 
tially distributed service time with the same mean, for the above three PH distributions we 
take three corresponding exponential distributions with service rates = 2.7500, /i(2) = 
3.4118 and ij,{3) = 2.3529, respectively. Table [9] illustrates how the doubly exponential 
solution (vTi to vTs) depends on the three service rates. Since the exponential distribution 
has a lower variance than the PH distribution, it is seen from Tables [8] and [9] that the 
service time has lower variance, 7r/c(Exp)< 7rfc(PH)e. 

Table 9: The doubly exponential solution depends on exponential service rates 





At(l) = 2.7500 


Ai(2) = 3.4118 


m(3) = 2.3529 




0.3636 


0.2931 


0.4250 


7r2 


0.0481 


0.0252 


0.0768 




8.408e-04 


1.858C-04 


0.0025 


7r4 


2.571e-07 


1.012C-08 


2.667C-06 




2.402e-14 


3.004C-17 


3.030e-12 



For the PH and exponential service times, the following two figures provides a compar- 
ison for the expected sojourn time. Clearly, the PH service time makes the lower expected 
sojourn time. 
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For m = 3,£i = 5,a(l) = (1/3,1/3,1/3) and a (2) = (1/12,7/12,1/3), 

/ -10 2 4 \ 

T= 3 _7 4 

\ 2 -5 / 

Table [10] shows how the doubly exponential solution (vri to iTi) depends on the vectors 
a (1) and a (2), respectively. 



Table 10: The doubly exponential solution depends on the vectors a 







12 ' 12 ' 3 ' 




(0.0741, 0.1358 , 0.2346) 


(0.0602, 0.1728, 0.2531) 


"■2 


(5.619e-05, 1.030e-05, 1.779e-04 ) 


(7.182e-05, 2.063e-04, 3.020e-04) 


1T3 


(1.411e-20, 2.587e-20, 4.469c-20) 


(1.739C-19, 4.993C-19, 7.311c-19) 


1T4 


(1.410e-98, 2.586e-98, 4.466e-98) 


(1.444e-92, 4.148e-92, 6.074e-92) 



6 Concluding remarks 

In this paper, we provide a matrix-analytic solution for supermarket models. We describe 
the supermarket model with PH service times as a system of differential vector equations, 
and provide a doubly exponential solution to the fixed point of the system of differential 
vector equations. We also provide some numerical examples to illustrate that our ap- 
proach is effective and efficient in the study of randomized load balancing schemes with 
non-exponential service requirements, such as, Erlang service time distributions, hyper- 
exponential service time distributions and PH service time distributions. We expect that 
this approach will be applicable to study other randomized load balancing schemes, for 
example, generalizing the arrival process to non-Poisson such as renewal process or Marko- 
vian arrival process, generalizing the service times to general probability distributions, and 
analyzing retrial and processor-sharing service disciplines. 
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